Load required packages for entire analysis
dlnm.pkg <- c("tidyverse", "lubridate","xts","ggthemes","PerformanceAnalytics","reshape2","rugarch","timetk","parallel","timeSeries","tseries","data.table","ggplot2","dlnm","broom","caret","gridExtra","splines","splines2","pspline","cowplot","mgcv","spi","chron","gridGraphics","grid","pscl","MASS", "AER", "Hmisc", "MuMIn", "VGAM", "forecast", "seasonal", "plotly", "ggmap", "rgeos", "tmap", "maptools", "maps", "ggfortify", "htmltools","webshot","knitr","flexdashboard", "imager", "httr", "curl", "gmodels","jpeg","rgdal","raster","sp","here")

lapply(dlnm.pkg, library, character.only=TRUE)
Load Malawi map dataset
dlnmtmp <- tempfile()
download.file("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/malawi_map.zip", destfile=dlnmtmp)
unzip(dlnmtmp, exdir = ".")
malawi.map <- rgdal::readOGR(".","malawi_map")
Subset Blantyre map from Malawi map
blantyre1.map <- malawi.map@data$OBJECTID >289 & malawi.map@data$OBJECTID <297 
blantyre2.map <- malawi.map@data$OBJECTID >308 & malawi.map@data$OBJECTID <311
blantyre3.map <- malawi.map@data$OBJECTID >342
Convert shape files into dataframe
blantyre.map <- rbind(fortify(malawi.map[blantyre1.map,]), fortify(malawi.map[blantyre2.map,]), fortify(malawi.map[blantyre3.map,]))
blantyre.map$id <- as.integer(blantyre.map$id)
Load Blantyre demographics’ and features’ dataset
blantyre.demog <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_demog.csv"))
    map.features <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_features.csv"))
    blantyre.demog$id <- as.integer(blantyre.demog$id)
    map.all <- merge(x=blantyre.map, y=blantyre.demog, by="id", x.all=TRUE)
    rm(list = ls()[grep("^blantyre", ls())])
Plot the map of Blantyre (Figure 1)
ggplot() + 
geom_polygon(data=map.all,aes(x=long,y=lat,group=group,fill=popc),colour="gray50") + 
theme_classic() + 
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) + 
labs(fill="(1998 - 2008) Population censuses") + 
xlab("Longitude") + 
ylab("Latitude") + 
geom_point(data=map.features, aes(x=long, y=lat, shape=Geolocation, size=Geolocation), color="black") +
scale_shape_manual(values=c(17, 16, 3)) +
scale_size_manual(values=c(2,4,3)) + 
theme(legend.key.height=unit(0.8, "line")) + 
theme(legend.key.width=unit(0.8, "line"))

Load case dataset and subset by serovar
case <-read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/case.csv"))
case$case_date <- dmy(case$case_date)
case$case_count <- c(1)
case.typhi <- subset(case, organism=="typhi")
case.iNTS <- subset(case, organism=="iNTS")
Assign 0 when no case was diagnosed on that day
case.typhi <-aggregate(case.typhi$case_count, by=list(case.typhi$case_date), FUN=sum, na.rm=TRUE)
names(case.typhi) <- c("date", "case_count")
case.typhi <- merge(case.typhi, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.typhi[is.na(case.typhi)]<-0

case.iNTS <-aggregate(case.iNTS$case_count, by=list(case.iNTS$case_date), FUN=sum, na.rm=TRUE)
names(case.iNTS) <- c("date", "case_count")
case.iNTS <- merge(case.iNTS, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.iNTS[is.na(case.iNTS)]<-0
Convert to XTS and aggregate monthly
case.typhi = as.xts(case.typhi[,-1,drop=FALSE], order.by=as.Date(case.typhi[,1]))
case.iNTS = as.xts(case.iNTS[,-1,drop=FALSE], order.by=as.Date(case.iNTS[,1]))
case.typhi <- apply.monthly(case.typhi, FUN=sum)
case.iNTS <- apply.monthly(case.iNTS, FUN=sum)
Load Blantyre climate dataset
climate <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/climate.csv"))
Define average rain/temp between chileka & chichiri stations
climate$climate_date <- dmy(climate$date)
climate$rainfall <- (climate$chil_r+climate$chic_r)/2 
climate$temperature <- ((climate$chil_mint+climate$chil_maxt)/2+(climate$chic_mint+climate$chic_maxt)/2)/2
climate <- subset(climate, select=c(climate_date, rainfall, temperature))
climate.rain <- subset(climate, select=c(climate_date, rainfall))
climate.temp <- subset(climate, select=c(climate_date, temperature))
Convert to XTS and aggregate monthly
climate.rain <- as.xts(climate.rain[,-1, drop=FALSE], order.by=as.Date(climate.rain[,1]))
climate.temp <- as.xts(climate.temp[,-1, drop=FALSE], order.by=as.Date(climate.temp[,1]))
climate.rain <- apply.monthly(climate.rain, FUN=mean)
climate.temp <- apply.monthly(climate.temp, FUN=mean)
Create datasets for time series decomposition
case.typhi <- tk_tbl(case.typhi, preserve_index=TRUE, rename_index="date") 
case.iNTS <- tk_tbl(case.iNTS, preserve_index=TRUE, rename_index="date") 
climate.rain <- tk_tbl(climate.rain, preserve_index=TRUE, rename_index="date") 
climate.temp <- tk_tbl(climate.temp, preserve_index=TRUE, rename_index="date")
Seasonal-adjusted NTS cases
case.iNTS.ts <- ts(na.omit(case.iNTS$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.iNTS.ts)))), preserve_index=FALSE)
case.iNTS$case_count_sea <- trend_n$value 
case.iNTS$case_count_sea[case.iNTS$case_count_sea < 0] <- 0
setnames(case.iNTS, old="case_count", new="case_count_obs")
trend_n <-tk_tbl(exp(mstl(log(case.iNTS.ts))))
case.iNTS$case_count_tre <- trend_n$Trend
Seasonal-adjusted typhoid cases
case.typhi.ts <- ts(na.omit(case.typhi$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.typhi.ts+1)))), preserve_index=FALSE)
case.typhi$case_count_sea <-trend_n$value
case.typhi$case_count_sea[case.typhi$case_count_sea < 0] <- 0
setnames(case.typhi, old="case_count", new="case_count_obs")
trend_n <- tk_tbl(exp(mstl(log(case.typhi.ts+1))))
case.typhi$case_count_tre <- trend_n$Trend
Seasonal-adjusted rainfall
climate.rain.ts <- ts(na.omit(climate.rain$rainfall), frequency=12)
trend_n <- tk_tbl(seasadj(mstl(climate.rain.ts)), preserve_index=FALSE)
climate.rain$rainfall_sea <- trend_n$value
climate.rain$rainfall_sea[climate.rain$rainfall_sea < 0] <- 0
setnames(climate.rain, old="rainfall", new="rainfall_obs")
trend_n <-tk_tbl(mstl(climate.rain.ts))
climate.rain$rain_tre <- trend_n$Trend
Seasonal-adjusted temperature
climate.temp.ts <- ts(na.omit(climate.temp$temperature), frequency=12)
trend_n <-tk_tbl(seasadj(mstl(climate.temp.ts)), preserve_index=FALSE)
climate.temp$temperature_sea <- trend_n$value
climate.temp$temperature_sea[climate.temp$temperature_sea < 0] <- 0
setnames(climate.temp, old="temperature", new="temperature_obs")
trend_n <- tk_tbl(mstl(climate.temp.ts))
climate.temp$temp_tre <- trend_n$Trend

rm(list=ls()[grep("^trend_n", ls())])
Plots of decomposed series (Supplementary Figure 1)
x<-case.iNTS.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="A",xlab="Years (2000-2015)",size=1,colour="orange2",is.date=FALSE)+theme_bw()
y<-case.typhi.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="B",xlab="Years (2000-2015)",size=1,colour="red2",is.date=FALSE)+theme_bw()
z<-climate.rain.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="C",xlab="Years (2000-2015)",size=1,colour="blue2",is.date=FALSE)+theme_bw()
v<-climate.temp.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="D",xlab="Years (2000-2015)",size=1,colour="green2",is.date=FALSE)+theme_bw()
grid.arrange(grobs=list(x,y,z,v), ncol=4, nrow=1)

Plots of seasonal-adjusted incidences (Supplementary Figure 2)
Th <- theme(plot.title=element_text(hjust=0)) + 
theme(axis.title.x=element_text(size=10)) + 
theme(axis.title.y=element_text(size=10)) +
theme(axis.text.x=element_text(face="bold", size=10), axis.text.y=element_text(face="bold", size=10)) + 
theme(legend.key.height=unit(1, "line")) + 
theme(legend.key.width=unit(1, "line")) 

S1 <- ggplot(as.data.frame(case.iNTS)) + 
geom_line(aes(date, case_count_obs, color="Original"), size=0.8) + 
geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) + 
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="orange2")) +
labs(title="A", x="Year", y="NTS cases") + Th +
theme(legend.justification=c(0.5,0), legend.position=c(0.7,0.7), legend.text=element_text(size=10), 
legend.title=element_text(size=0))

S2<-ggplot(as.data.frame(case.typhi)) + 
geom_line(aes(date, case_count_obs, color="Original"), size=0.8) + 
geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) + 
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="red2")) +
labs(title="B", x ="Year", y="Typhoid cases") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10), 
legend.title = element_text(size=0))

S3<-ggplot(as.data.frame(climate.rain)) + 
geom_line(aes(date, rainfall_obs, color="Original"), size=0.8) + 
geom_line(aes(date, rainfall_sea, color="Seasonal-adjusted"), size=0.8) + 
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="blue2")) +
labs(title="C", x ="Year", y = "Rainfall (mm)") + Th + 
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10), 
legend.title = element_text(size=0))

S4<-ggplot(as.data.frame(climate.temp)) + 
geom_line(aes(date, temperature_obs, color="Original"), size=0.8) + 
geom_line(aes(date, temperature_sea, color="Seasonal-adjusted"), size=0.8) + 
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="green2")) + 
labs(title="D", x ="Year", y = "Temperature (°C)") + Th + 
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10), 
legend.title = element_text(size=0))

grid.arrange(grobs=list(S1, S2), ncol=2, nrow=1)

Inter/extraporate population denominators
census.year <-c(1998, 2008)
census.popn <-c(809397, 1022680)
census.count.1998.2008 <-approx(census.year, census.popn, n=11) 
census.count.2009.2015 <-approxExtrap(census.year, census.popn, xout=c(2009, 2010, 2011, 2012, 2013, 2014, 2015))
Calculate NTS/typhoid incidence rates
case.iNTS$census[year(case.iNTS$date)==2000]<-852054; case.iNTS$census[year(case.iNTS$date)==2001]<-873382
case.iNTS$census[year(case.iNTS$date)==2002]<-894710; case.iNTS$census[year(case.iNTS$date)==2003]<-916039
case.iNTS$census[year(case.iNTS$date)==2004]<-937367; case.iNTS$census[year(case.iNTS$date)==2005]<-958695
case.iNTS$census[year(case.iNTS$date)==2006]<-980023; case.iNTS$census[year(case.iNTS$date)==2007]<-1001352
case.iNTS$census[year(case.iNTS$date)==2008]<-1022680; case.iNTS$census[year(case.iNTS$date)==2009]<-1044008
case.iNTS$census[year(case.iNTS$date)==2010]<-1065337; case.iNTS$census[year(case.iNTS$date)==2011]<-1086665
case.iNTS$census[year(case.iNTS$date)==2012]<-1107993; case.iNTS$census[year(case.iNTS$date)==2013]<-1129322
case.iNTS$census[year(case.iNTS$date)==2014]<-1150650; case.iNTS$census[year(case.iNTS$date)==2015]<-1171978
case.iNTS$incid_sea <-case.iNTS$case_count_sea*100000/case.iNTS$census 
case.iNTS$incid_obs <-case.iNTS$case_count_obs*100000/case.iNTS$census 

case.typhi$census[year(case.typhi$date)==2000]<-852054; case.typhi$census[year(case.typhi$date)==2001]<-873382
case.typhi$census[year(case.typhi$date)==2002]<-894710; case.typhi$census[year(case.typhi$date)==2003]<-916039
case.typhi$census[year(case.typhi$date)==2004]<-937367; case.typhi$census[year(case.typhi$date)==2005]<-958695
case.typhi$census[year(case.typhi$date)==2006]<-980023; case.typhi$census[year(case.typhi$date)==2007]<-1001352
case.typhi$census[year(case.typhi$date)==2008]<-1022680; case.typhi$census[year(case.typhi$date)==2009]<-1044008
case.typhi$census[year(case.typhi$date)==2010]<-1065337; case.typhi$census[year(case.typhi$date)==2011]<-1086665
case.typhi$census[year(case.typhi$date)==2012]<-1107993; case.typhi$census[year(case.typhi$date)==2013]<-1129322
case.typhi$census[year(case.typhi$date)==2014]<-1150650; case.typhi$census[year(case.typhi$date)==2015]<-1171978
case.typhi$incid_sea <-case.typhi$case_count_sea*100000/case.typhi$census 
case.typhi$incid_obs <-case.typhi$case_count_obs*100000/case.typhi$census
Plots of Age-sex distribution of cases (Figure 3)
case$sex[case$sex == ""] <- NA
case$age[case$age == ""] <- NA
case$date <- ymd(case$case_date)
case$year <- year(case$date)

agesex.p1 <- ggplot(subset(case, organism=="iNTS" & !is.na(age) & sex != "Unknown"), aes(x=age, color=sex)) + 
geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) + 
scale_color_manual(values=c(Male="gray40",Female="red3")) + 
theme_bw() + 
scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) + 
labs(title="A", x="Age (years)", y="Number of NTS cases") + 
theme(plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) + 
theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) + 
labs(color="Missing sex: 2,596 (32.3%)") + 
theme(legend.key.height=unit(1,"line")) + 
theme(legend.key.width=unit(1,"line"))

agesex.p2<- ggplot(subset(case, organism=="typhi" & !is.na(age) & sex !="Unknown"), aes(x=age, color=sex)) + 
geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) + 
scale_color_manual(values=c(Male="gray40",Female="red3")) + 
theme_bw() + 
scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) + 
labs(title="B", x="Age (years)", y="Number of typhoid cases") + 
theme(plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) + 
theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) + 
labs(color="Missing sex: 61 (2.4%)") + 
theme(legend.key.height=unit(1,"line")) + 
theme(legend.key.width=unit(1,"line"))

grid.arrange(grobs=list(agesex.p1, agesex.p2), ncol=2, nrow=1)

Generating yearly and monthly dynamics
case.iNTS.spi <- subset(case.iNTS, year(case.iNTS$date)<2011, select=c(date,incid_sea)) 
case.iNTS.spi$month <-month(case.iNTS.spi$date)
case.iNTS.spi$year <- year(case.iNTS.spi$date); case.iNTS.spi$date <- NULL
case.iNTS.spi <- spread(case.iNTS.spi, year, incid_sea); case.iNTS.spi <- as.matrix(case.iNTS.spi)[,-1]
YMD1 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~case.iNTS.spi, type="contour", colorscale='jet', contours=list(showlabels=TRUE)) %>% 
colorbar(title="NTS incidence per \n 100,000 population") %>%
layout(title="A", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
case.typhi.spi <- subset(case.typhi, year(case.typhi$date)>2010, select=c(date,incid_sea)) 
case.typhi.spi$month <- month(case.typhi.spi$date)
case.typhi.spi$year <- year(case.typhi.spi$date); case.typhi.spi$date <- NULL
case.typhi.spi <- spread(case.typhi.spi, year, incid_sea); case.typhi.spi <- as.matrix(case.typhi.spi)[,-1]
YMD2 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~case.typhi.spi, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% 
colorbar(title="Typhoid incidence per \n 100,000 population") %>%
layout(title="D", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spin <- subset(climate.rain, year(climate.rain$date)<2011, select=c(date,rainfall_obs)) 
climate.rain.spin$month <- month(climate.rain.spin$date)
climate.rain.spin$year <- year(climate.rain.spin$date); climate.rain.spin$date <- NULL
climate.rain.spin <- spread(climate.rain.spin, year, rainfall_obs); climate.rain.spin <- as.matrix(climate.rain.spin)[,-1]
YMD3 <- plot_ly(x = c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.rain.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% colorbar(title="Rainfall (mm)") %>% 
layout(title="B", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spit <- subset(climate.rain, year(climate.rain$date)>2010, select=c(date,rainfall_obs)) 
climate.rain.spit$month <- month(climate.rain.spit$date)
climate.rain.spit$year <- year(climate.rain.spit$date); climate.rain.spit$date <- NULL
climate.rain.spit <- spread(climate.rain.spit, year, rainfall_obs); climate.rain.spit <- as.matrix(climate.rain.spit)[,-1]
YMD4 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.rain.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% 
colorbar(title="Rainfall (mm)") %>%
layout(title="E", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spin <- subset(climate.temp, year(climate.temp$date)<2011, select=c(date,temperature_obs)) 
climate.temp.spin$month <- month(climate.temp.spin$date)
climate.temp.spin$year <- year(climate.temp.spin$date); climate.temp.spin$date <- NULL
climate.temp.spin <- spread(climate.temp.spin, year, temperature_obs); climate.temp.spin <- as.matrix(climate.temp.spin)[,-1]
YMD5 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.temp.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% 
colorbar(title="Temperature (°C)") %>%
layout(title="C", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spit <- subset(climate.temp, year(climate.temp$date)>2010, select=c(date,temperature_obs)) 
climate.temp.spit$month <- month(climate.temp.spit$date)
climate.temp.spit$year <- year(climate.temp.spit$date); climate.temp.spit$date <- NULL
climate.temp.spit <- spread(climate.temp.spit, year, temperature_obs); climate.temp.spit <- as.matrix(climate.temp.spit)[,-1]
YMD6 <- plot_ly(x=c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.temp.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Temperature (°C)") %>%
layout(title="F", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size = 13))
Plots of yearly and monthly dynamics (Figure 2)
htmltools::tags$div(
style = "display: flex; flex-wrap: wrap",
tags$div(YMD1, style = "width: 33%;"),
tags$div(YMD3, style = "width: 33%;"),
tags$div(YMD5, style = "width: 33%;"),
tags$div(YMD2, style = "width: 33%;"),
tags$div(YMD4, style = "width: 33%;"),
tags$div(YMD6, style = "width: 33%;")
)
Constructing DLNM datasets for NTS and Typhoid
mo.dlnmN <- bind_cols(case.iNTS, climate.rain, climate.temp, id=NULL)
mo.dlnmN$date1 <- mo.dlnmN$date2 <- NULL
mo.dlnmN <- subset(mo.dlnmN, year(date) < 2011)
mo.dlnmN$time <- seq.int(from = 1, to=132, by=1)
mo.dlnmN$year <- year(mo.dlnmN$date)
mo.dlnmN$month <- month(mo.dlnmN$date)
mo.dlnmN$incid_seaX<-round(mo.dlnmN$incid_sea, digits = 0)
mo.dlnmN$incid_obsX<-round(mo.dlnmN$incid_obs, digits = 0)

mo.dlnmT <- bind_cols(case.typhi, climate.rain, climate.temp, id=NULL)
mo.dlnmT$date1 <- mo.dlnmT$date2 <- NULL
mo.dlnmT <- subset(mo.dlnmT, year(date) > 2010)
mo.dlnmT$time <- seq.int(from = 1, to=60, by=1)
mo.dlnmT$year <- year(mo.dlnmT$date)
mo.dlnmT$month <- month(mo.dlnmT$date)
mo.dlnmT$incid_seaX<-round(mo.dlnmT$incid_sea, digits = 0)
Switching from AIC to QAIC for NTS and Typhoid
nts.quasipoisson <- function(...) { 
  res <- quasipoisson(...)
  res$aic <- poisson(...)$aic 
  res
}
typ.quasipoisson <- function(...) { 
  res <- quasipoisson(...) 
  res$aic <- poisson(...)$aic 
  res
}
Calculate QAIC values to select optimal models (Supplementary Table)
QAICtable <- data.frame(model.no=rep(NA,27), lag.df=rep(NA,27), fx1.df=rep(NA,27), fx2.df=rep(NA,27), 
                    QAIC.ntsR=rep(NA,27), QAIC.ntsT=rep(NA,27), QAIC.nts=rep(NA,27), QAIC.typR=rep(NA,27), 
                    QAIC.typT=rep(NA,27), QAIC.typ=rep(NA,27))
l=1
for(k in 3:5){
for(j in 3:5){
  for(i in 3:5){
nts.lagknots <- logknots(8, fun="ns", df=i)
nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=j)
nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=k)
nts.mo.cb.rain <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots))
nts.mo.cb.temp <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots))
nts.modelR <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
nts.modelT <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.temp + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
nts.model <-  glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + nts.mo.cb.temp + year, family = nts.quasipoisson(), na.action=na.delete, mo.dlnmN)

typ.lagknots <- logknots(8, fun="ns", df=i)
typ.varknots.r=equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=j)
typ.varknots.t=equalknots(mo.dlnmT$temperature_obs, fun="ns", df=k)
typ.mo.cb.rain <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r), arglag=list(knots=typ.lagknots))
typ.mo.cb.temp <- crossbasis(mo.dlnmT$temperature_obs, lag =8, argvar = list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.model <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)

QAICtable[l,] <- c(l,i,j,k, QAIC(nts.modelR, chat=summary(nts.modelR)$dispersion), QAIC(nts.modelT, chat=summary(nts.modelT)$dispersion), 
                   QAIC(nts.model, chat=summary(nts.model)$dispersion), QAIC(typ.modelR, chat=summary(typ.modelR)$dispersion), QAIC(typ.modelT, chat=summary(typ.modelT)$dispersion), 
                   QAIC(typ.model, chat=1))
l=l+1  
  }
}
  }
kable(QAICtable)
model.no lag.df fx1.df fx2.df QAIC.ntsR QAIC.ntsT QAIC.nts QAIC.typR QAIC.typT QAIC.typ
1 3 3 3 570.8194 565.8631 562.3950 218.5004 220.2006 212.6658
2 4 3 3 576.3831 571.4838 573.3630 223.1714 216.1853 220.6390
3 5 3 3 581.4852 574.9872 583.4617 228.2399 221.1698 230.8468
4 3 4 3 571.1280 565.8631 565.3731 211.3436 220.2006 216.2219
5 4 4 3 578.3398 571.4838 578.1025 218.1843 216.1853 228.0142
6 5 4 3 585.4681 574.9872 590.3358 222.7394 221.1698 237.6362
7 3 5 3 572.7386 565.8631 566.1964 216.3176 220.2006 221.5570
8 4 5 3 581.6967 571.4838 581.3661 222.4327 216.1853 234.0871
9 5 5 3 590.5606 574.9872 594.8792 227.8628 221.1698 246.2023
10 3 3 4 570.8194 571.7243 566.2250 218.5004 220.7893 215.5219
11 4 3 4 576.3831 577.5401 578.9833 223.1714 220.2852 225.9024
12 5 3 4 581.4852 583.4447 589.9669 228.2399 226.2618 236.0380
13 3 4 4 571.1280 571.7243 568.6740 211.3436 220.7893 220.0702
14 4 4 4 578.3398 577.5401 583.3435 218.1843 220.2852 233.1946
15 5 4 4 585.4681 583.4447 597.5272 222.7394 226.2618 245.6773
16 3 5 4 572.7386 571.7243 571.2271 216.3176 220.7893 224.5520
17 4 5 4 581.6967 577.5401 588.2408 222.4327 220.2852 239.9495
18 5 5 4 590.5606 583.4447 603.1850 227.8628 226.2618 252.6551
19 3 3 5 570.8194 574.0369 566.7772 218.5004 220.7423 220.4298
20 4 3 5 576.3831 580.7484 580.7418 223.1714 220.9285 232.7209
21 5 3 5 581.4852 587.7540 593.8723 228.2399 225.5709 245.1746
22 3 4 5 571.1280 574.0369 569.9469 211.3436 220.7423 224.8192
23 4 4 5 578.3398 580.7484 586.2703 218.1843 220.9285 239.8864
24 5 4 5 585.4681 587.7540 601.6447 222.7394 225.5709 254.6679
25 3 5 5 572.7386 574.0369 575.0626 216.3176 220.7423 229.8974
26 4 5 5 581.6967 580.7484 593.7784 222.4327 220.9285 246.7405
27 5 5 5 590.5606 587.7540 610.6531 227.8628 225.5709 261.6885
Reformulate cross-basis matrix using selected best NTS model
sort(mo.dlnmN$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmN$temperature_obs, decreasing=FALSE)
nts.lagknots.r <- logknots(8, fun="ns", df=3)
nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=3)
nts.lagknots.t <- logknots(8, fun="ns", df=3)
nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=3)
nts.cb.rain <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots.r))
nts.cb.temp <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots.t))
summary(nts.cb.rain)
summary(nts.cb.temp)
Fit NTS model
nts.model <- glm(mo.dlnmN$incid_seaX ~ nts.cb.rain + nts.cb.temp + year, family=quasipoisson(), na.action=na.delete, mo.dlnmN)
Investigate NTS deviance residuals
par(mfrow=c(1,2))
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Original model with Rain & Temp", xlim=c(0,8))
nts.model <- update(nts.model,.~.+Lag(residuals(nts.model,type="deviance"), 1)) #add residuals at lag 1 
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Adjusted model with Rain & Temp", xlim=c(0,8))

Run NTS predictions
nts.pred.rain <- crosspred(nts.cb.rain, nts.model, cen=0, by=0.2)
nts.pred.temp <- crosspred(nts.cb.temp, nts.model, cen=23, by=0.2)
Plots of NTS predictions due to rainfall (Figure 4)
#plot(nts.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1.1, cex.axis=1.1, col="gray80",main="A")
plot(nts.pred.rain, "contour", key.title=title("NTS"), plot.title=title("", xlab ="Rainfall (mm)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="A"))

plot(nts.pred.rain, "slices", xlab="Lag (month) [given 9 mm] ", var=c(9), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="B")

plot(nts.pred.rain, "slices", xlab="Lag (month) [given 13 mm] ", var=c(13), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="C")

Plots of NTS predictions due to temperature (Figure 6)
#plot(nts.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1, cex.axis=1, col="gray80",main="A")
plot(nts.pred.temp, "contour", key.title=title("NTS"), plot.title=title("", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="A"))

plot(nts.pred.temp, xlab="Lag (month) [given 19 °C]", "slices",var=c(19), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="B")

plot(nts.pred.temp, xlab="Lag (month) [given 29 °C]", "slices",var=c(29), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="C")

Reformulate cross-basis matrix using selected best Typhoid models
sort(mo.dlnmT$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmT$temperature_obs, decreasing=FALSE)
typ.lagknots <- logknots(8, fun="ns", df=3)
typ.varknots.r1 <- equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=4)
typ.varknots.r2 <- equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=3)
typ.varknots.t <- equalknots(mo.dlnmT$temperature_obs, fun="ns", df=3)
typ.cb.rain1 <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r1), arglag=list(knots=typ.lagknots))
typ.cb.rain2 <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r2), arglag=list(knots=typ.lagknots))
typ.cb.temp <- crossbasis(mo.dlnmT$temperature_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
summary(typ.cb.rain1)
summary(typ.cb.rain2)
summary(typ.cb.temp)
Fit Typhoid models
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.cb.rain1 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.cb.temp + typ.cb.rain2 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
Investigate deviance residuals for fitted tyhpoid models
par(mfrow=c(1,3))
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="Original model with Rain",xlim=c(0,8))
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Original model with Rain & Temp",xlim=c(0,8))
typ.modelT <- update(typ.modelT,.~.+Lag(residuals(typ.modelT,type="deviance"),1)) #add residuals at lag 1
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Adjusted model with Rain & Temp",xlim=c(0,8))

Run typhoid predictions
typ.pred.rain <- crosspred(typ.cb.rain1, typ.modelR, cen=0, by=0.2)
typ.pred.temp <- crosspred(typ.cb.temp, typ.modelT, cen=23, by=0.2)
Plots of TYP predictions due to rainfall (Figure 5)
#plot(typ.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of typhoid", theta=210, phi=10, lphi=150, col="gray80", main="D", cex.lab=1.1, cex.axis=1.1)
plot(typ.pred.rain, "contour", key.title=title("TYP"), plot.title=title("", xlab ="Rainfall (mm)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="D"))

plot(typ.pred.rain, xlab="Lag (month) [given 9 mm]", "slices", var=c(9), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="E")

plot(typ.pred.rain, xlab="Lag (month) [given 16 mm]", "slices", var=c(16), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="F")

Plots of TYP predictions due to temperature (Figure 7)
#plot(typ.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of typhoid", theta=150, phi=5, lphi=100, col="gray80", main="D", cex.lab=1.1, cex.axis=1.1)
plot(typ.pred.temp, "contour", key.title=title("TYP"), plot.title=title("", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="D"))

plot(typ.pred.temp, xlab="Lag (month) [given 19 °C]", "slices", var=c(19), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="E")

plot(typ.pred.temp, xlab="Lag (month) [given 25 °C]", "slices", var=c(25), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="F")

Plots of validation checks (Supplementary Figure S4)
par(mfrow=c(3,4))
plot(mo.dlnmN$date, residuals(nts.model,type="deviance"), pch=19, cex=0.8, col=grey(0.4),main="A", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2) 
abline(h=0,lty=2,lwd=3)
pacf(residuals(nts.model,type="deviance"),na.action=na.omit,main="B",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
acf(residuals(nts.model,type="deviance"),na.action=na.omit,main="C",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
plot(matrix(mo.dlnmN$incid_seaX),matrix(predict(nts.model, type="response")), col=c("orange2","gray2"),main="D",xlab="Predicted NTS incidence", ylab="Observed NTS incidence",pch=19,cex=0.8,cex.lab=1.2,cex.axis=1.2)
legend("topleft", legend=c("observed", "predicted"), col=c("grey30", "orange2"), cex=0.8, pch=19)

plot(mo.dlnmT$date, residuals(typ.modelR,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="E", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="F",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="G",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
plot(matrix(mo.dlnmT$incid_seaX),matrix(predict(typ.modelR, type="response")), col=c("red2","gray2"),main="H",xlab="Predicted typhoid incidence", ylab="Observed typhoid incidence",pch=19,cex=0.8,cex.lab=1.2,cex.axis=1.2)
legend("topleft", legend=c("observed", "predicted"), col=c("grey30", "red2"), cex=0.8, pch=19)

plot(mo.dlnmT$date, residuals(typ.modelT,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="I", ylab="Residuals", xlab="Year",cex.lab=1.2,cex.axis=1.2) 
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="J",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="K",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
plot(matrix(mo.dlnmT$incid_seaX),matrix(predict(typ.modelT, type="response")), col=c("red2","gray2"),main="L",xlab="Predicted typhoid incidence", ylab="Observed typhoid incidence",pch=19,cex=0.8,cex.lab=1.2,cex.axis=1.2)
legend("topleft", legend=c("observed", "predicted"), col=c("grey30", "red2"), cex=0.8, pch=19)

Plots of predictions due to rainfall/temperature (Table 1)
RRTable <- data.frame(rbind(nts.pred.rain$matRRfit["9",], nts.pred.rain$matRRlow["9",], nts.pred.rain$matRRhigh["9",],
                 nts.pred.rain$matRRfit["13",],nts.pred.rain$matRRlow["13",],nts.pred.rain$matRRhigh["13",],
                 nts.pred.temp$matRRfit["19",],nts.pred.temp$matRRlow["19",],nts.pred.temp$matRRhigh["19",],
                 nts.pred.temp$matRRfit["29",],nts.pred.temp$matRRlow["29",],nts.pred.temp$matRRhigh["29",],
                 typ.pred.rain$matRRfit["9",],typ.pred.rain$matRRlow["9",],typ.pred.rain$matRRhigh["9",],
                 typ.pred.rain$matRRfit["16",],typ.pred.rain$matRRlow["16",],typ.pred.rain$matRRhigh["16",],
                 typ.pred.temp$matRRfit["19",],typ.pred.temp$matRRlow["19",],typ.pred.temp$matRRhigh["19",],
                 typ.pred.temp$matRRfit["25",],typ.pred.temp$matRRlow["25",],typ.pred.temp$matRRhigh["25",]
                 ))
kable(RRTable)
lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7 lag8
1.2248467 1.2514626 1.2554780 1.2238951 1.1635420 1.0841819 0.9951394 0.9042881 0.8176177
1.0254867 1.1162267 1.1453653 1.1143857 1.0596817 0.9906237 0.9027750 0.7975713 0.6897005
1.4629635 1.4030831 1.3761766 1.3441658 1.2775817 1.1865762 1.0969538 1.0252839 0.9692593
1.1817325 1.0655402 0.9723354 0.9040825 0.8545278 0.8183618 0.7914806 0.7705230 0.7525862
0.9099591 0.8927064 0.8190710 0.7501790 0.7124350 0.6968964 0.6814001 0.6417570 0.5822991
1.5346752 1.2718359 1.1542787 1.0895601 1.0249607 0.9609980 0.9193445 0.9251255 0.9726719
1.0822001 1.1309828 1.1629884 1.1659357 1.1432279 1.1012306 1.0467408 0.9861515 0.9249535
0.9602616 1.0318836 1.0674043 1.0690932 1.0514755 1.0192753 0.9708904 0.9049599 0.8287981
1.2196229 1.2395991 1.2671318 1.2715505 1.2429867 1.1897756 1.1285169 1.0746274 1.0322647
0.7825270 0.7912493 0.7912237 0.7775272 0.7525056 0.7194590 0.6815991 0.6418050 0.6024941
0.6253076 0.6887091 0.6928915 0.6692501 0.6487554 0.6292260 0.5946514 0.5350856 0.4650197
0.9792756 0.9090565 0.9035107 0.9033223 0.8728476 0.8226317 0.7812600 0.7698090 0.7806102
0.9870340 1.2979004 1.6063110 1.8078055 1.8722857 1.8143089 1.6725894 1.4915122 1.3081090
0.6284634 0.9438326 1.2552773 1.4434731 1.5069006 1.4634299 1.3249680 1.1234975 0.9124989
1.5501876 1.7847926 2.0555100 2.2640954 2.3262675 2.2493164 2.1114137 1.9800744 1.8752341
0.2199818 0.2071173 0.2101761 0.2398628 0.3033687 0.4165633 0.6083705 0.9257774 1.4380388
0.0679390 0.0908122 0.1004582 0.1120403 0.1444326 0.2128310 0.3269832 0.4616799 0.5759679
0.7122860 0.4723769 0.4397253 0.5135131 0.6372006 0.8153182 1.1319073 1.8564025 3.5904008
1.4960476 1.5619360 1.5874626 1.5468052 1.4526163 1.3245023 1.1812571 1.0380766 0.9055470
1.0424914 1.1457962 1.0968536 1.0279449 0.9761052 0.9252186 0.8373032 0.6896236 0.5229565
2.1469324 2.1292129 2.2975149 2.3275628 2.1617487 1.8960993 1.6665030 1.5625959 1.5680377
0.6881750 0.9207072 1.1640344 1.3467161 1.4416997 1.4504554 1.3928588 1.2966573 1.1885051
0.4844865 0.7123388 0.8891649 0.9975957 1.0683040 1.1081654 1.1070709 1.0400232 0.9069586
0.9774984 1.1900260 1.5238750 1.8180153 1.9456054 1.8984720 1.7524222 1.6166181 1.5574519
Vary degrees of freedom for sensitivity analyses of NTS/typhoid predictions
nts.cb.rain <- list()
nts.cb.temp <- list()
typ.cb.rain <- list()
typ.cb.temp <- list()
sort(mo.dlnmN$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmN$temperature_obs, decreasing=FALSE)
sort(mo.dlnmT$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmT$temperature_obs, decreasing=FALSE)
l=1
for(k in 3:5){
  for(j in 3:5){
    for(i in 3:5){
nts.lagknots <- logknots(8, fun="ns", df=i)
nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=j)
nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=k)
nts.cb.rain[[l]] <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots))
nts.cb.temp[[l]] <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots))

typ.lagknots <- logknots(8, fun="ns", df=i)
typ.varknots.r=equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=j)
typ.varknots.t=equalknots(mo.dlnmT$temperature_obs, fun="ns", df=k)
typ.cb.rain[[l]] <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r), arglag=list(knots=typ.lagknots))
typ.cb.temp[[l]] <- crossbasis(mo.dlnmT$temperature_obs, lag =8, argvar = list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
  
  l=l+1  
  }
  }
}
Formulate alternative sensitivity models and plot predictions of NTS (Supplementay Figure S5)
nts.cb.rain.s1 <- nts.cb.rain[[2]]; nts.cb.temp.s1 <- nts.cb.temp[[2]]
nts.cb.rain.s2 <- nts.cb.rain[[3]]; nts.cb.temp.s2 <- nts.cb.temp[[3]]
nts.cb.rain.s3 <- nts.cb.rain[[4]]; nts.cb.temp.s3 <- nts.cb.temp[[4]]
nts.cb.rain.s4 <- nts.cb.rain[[7]]; nts.cb.temp.s4 <- nts.cb.temp[[7]]
nts.cb.rain.s5 <- nts.cb.rain[[10]]; nts.cb.temp.s5 <- nts.cb.temp[[10]]
nts.cb.rain.s6 <- nts.cb.rain[[19]]; nts.cb.temp.s6 <- nts.cb.temp[[19]]

nts.model.s1 <- glm(incid_seaX~nts.cb.rain.s1 + nts.cb.temp.s1 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s2 <- glm(incid_seaX~nts.cb.rain.s2 + nts.cb.temp.s2 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s3 <- glm(incid_seaX~nts.cb.rain.s3 + nts.cb.temp.s3 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s4 <- glm(incid_seaX~nts.cb.rain.s4 + nts.cb.temp.s4 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s5 <- glm(incid_seaX~nts.cb.rain.s5 + nts.cb.temp.s5 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s6 <- glm(incid_seaX~nts.cb.rain.s6 + nts.cb.temp.s6 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)

nts.model.s1 <- update(nts.model.s1,.~.+Lag(residuals(nts.model.s1,type="deviance"),1)) #add residuals at lag 1
nts.model.s2 <- update(nts.model.s2,.~.+Lag(residuals(nts.model.s2,type="deviance"),1)) #to correct for partial
nts.model.s3 <- update(nts.model.s3,.~.+Lag(residuals(nts.model.s3,type="deviance"),1)) #autocorrelation
nts.model.s4 <- update(nts.model.s4,.~.+Lag(residuals(nts.model.s4,type="deviance"),1)) 
nts.model.s5 <- update(nts.model.s5,.~.+Lag(residuals(nts.model.s5,type="deviance"),1)) 
nts.model.s6 <- update(nts.model.s6,.~.+Lag(residuals(nts.model.s6,type="deviance"),1)) 

nts.pred.rain.s1 <- crosspred(nts.cb.rain.s1, nts.model.s1, cen=0, by=0.2)
nts.pred.rain.s2 <- crosspred(nts.cb.rain.s2, nts.model.s2, cen=0, by=0.2)
nts.pred.rain.s3 <- crosspred(nts.cb.rain.s3, nts.model.s3, cen=0, by=0.2)
nts.pred.rain.s4 <- crosspred(nts.cb.rain.s4, nts.model.s4, cen=0, by=0.2)
nts.pred.rain.s5 <- crosspred(nts.cb.rain.s5, nts.model.s5, cen=0, by=0.2)
nts.pred.rain.s6 <- crosspred(nts.cb.rain.s6, nts.model.s6, cen=0, by=0.2)

nts.pred.temp.s1 <- crosspred(nts.cb.temp.s1, nts.model.s1, cen=23, by=0.2)
nts.pred.temp.s2 <- crosspred(nts.cb.temp.s2, nts.model.s2, cen=23, by=0.2)
nts.pred.temp.s3 <- crosspred(nts.cb.temp.s3, nts.model.s3, cen=23, by=0.2)
nts.pred.temp.s4 <- crosspred(nts.cb.temp.s4, nts.model.s4, cen=23, by=0.2)
nts.pred.temp.s5 <- crosspred(nts.cb.temp.s5, nts.model.s5, cen=23, by=0.2)
nts.pred.temp.s6 <- crosspred(nts.cb.temp.s6, nts.model.s6, cen=23, by=0.2)

par(mfrow=c(2,6))
plot(nts.pred.rain.s1,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="A")
plot(nts.pred.rain.s2,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="C")
plot(nts.pred.rain.s3,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="E")
plot(nts.pred.rain.s4,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="G")
plot(nts.pred.rain.s5,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="I")
plot(nts.pred.rain.s6,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="K")

plot(nts.pred.temp.s1,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="B")
plot(nts.pred.temp.s2,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="D")
plot(nts.pred.temp.s3,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="F")
plot(nts.pred.temp.s4,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="H")
plot(nts.pred.temp.s5,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="J")
plot(nts.pred.temp.s6,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="L")

Formulate alternative sensitivity models and plot predictions of Typhoid (Supplementay Figure S6)
typ.cb.rain.s1 <- typ.cb.rain[[2]]; typ.cb.temp.s1 <- typ.cb.temp[[2]]
typ.cb.rain.s2 <- typ.cb.rain[[3]]; typ.cb.temp.s2 <- typ.cb.temp[[3]]
typ.cb.rain.s3 <- typ.cb.rain[[4]]; typ.cb.temp.s3 <- typ.cb.temp[[4]]
typ.cb.rain.s4 <- typ.cb.rain[[7]]; typ.cb.temp.s4 <- typ.cb.temp[[7]]
typ.cb.rain.s5 <- typ.cb.rain[[10]]; typ.cb.temp.s5 <- typ.cb.temp[[10]]
typ.cb.rain.s6 <- typ.cb.rain[[19]]; typ.cb.temp.s6 <- typ.cb.temp[[19]]

typ.modelR.s1 <- glm(incid_seaX~typ.cb.rain.s1 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s2 <- glm(incid_seaX~typ.cb.rain.s2 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s3 <- glm(incid_seaX~typ.cb.rain.s3 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s4 <- glm(incid_seaX~typ.cb.rain.s4 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s5 <- glm(incid_seaX~typ.cb.rain.s5 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s6 <- glm(incid_seaX~typ.cb.rain.s6 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s1 <- glm(incid_seaX~typ.cb.rain.s1 + typ.cb.temp.s1 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s2 <- glm(incid_seaX~typ.cb.rain.s2 + typ.cb.temp.s2 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s3 <- glm(incid_seaX~typ.cb.rain.s3 + typ.cb.temp.s3 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s4 <- glm(incid_seaX~typ.cb.rain.s4 + typ.cb.temp.s4 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s5 <- glm(incid_seaX~typ.cb.rain.s5 + typ.cb.temp.s5 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s6 <- glm(incid_seaX~typ.cb.rain.s6 + typ.cb.temp.s6 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
    
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),1))
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),2))
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),3))
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),4))
typ.modelR.s2 <- update(typ.modelR.s2,.~.+Lag(residuals(typ.modelR.s2,type="deviance"),1)) 
typ.modelR.s2 <- update(typ.modelR.s2,.~.+Lag(residuals(typ.modelR.s2,type="deviance"),2))
typ.modelR.s2 <- update(typ.modelR.s2,.~.+Lag(residuals(typ.modelR.s2,type="deviance"),3))
typ.modelR.s5 <- update(typ.modelR.s5,.~.+Lag(residuals(typ.modelR.s5,type="deviance"),1)) 
typ.modelR.s5 <- update(typ.modelR.s5,.~.+Lag(residuals(typ.modelR.s5,type="deviance"),2))
typ.modelR.s5 <- update(typ.modelR.s5,.~.+Lag(residuals(typ.modelR.s5,type="deviance"),3))
typ.modelR.s6 <- update(typ.modelR.s6,.~.+Lag(residuals(typ.modelR.s6,type="deviance"),1)) 
typ.modelR.s6 <- update(typ.modelR.s6,.~.+Lag(residuals(typ.modelR.s6,type="deviance"),2)) 
typ.modelR.s6 <- update(typ.modelR.s6,.~.+Lag(residuals(typ.modelR.s6,type="deviance"),3)) 

typ.modelT.s1 <- update(typ.modelT.s1,.~.+Lag(residuals(typ.modelT.s1,type="deviance"),1))
typ.modelT.s1 <- update(typ.modelT.s1,.~.+Lag(residuals(typ.modelT.s1,type="deviance"),3))
typ.modelT.s2 <- update(typ.modelT.s2,.~.+Lag(residuals(typ.modelT.s2,type="deviance"),1)) 
typ.modelT.s3 <- update(typ.modelT.s3,.~.+Lag(residuals(typ.modelT.s3,type="deviance"),1)) 
typ.modelT.s4 <- update(typ.modelT.s4,.~.+Lag(residuals(typ.modelT.s4,type="deviance"),1)) 
typ.modelT.s4 <- update(typ.modelT.s4,.~.+Lag(residuals(typ.modelT.s4,type="deviance"),2)) 
typ.modelT.s5 <- update(typ.modelT.s5,.~.+Lag(residuals(typ.modelT.s5,type="deviance"),1)) 
typ.modelT.s5 <- update(typ.modelT.s5,.~.+Lag(residuals(typ.modelT.s5,type="deviance"),4))
typ.modelT.s5 <- update(typ.modelT.s5,.~.+Lag(residuals(typ.modelT.s5,type="deviance"),3))
typ.modelT.s6 <- update(typ.modelT.s6,.~.+Lag(residuals(typ.modelT.s6,type="deviance"),1)) 

typ.pred.rain.s1 <- crosspred(typ.cb.rain.s1, typ.modelR.s1, cen=0, by=0.2)
typ.pred.rain.s2 <- crosspred(typ.cb.rain.s2, typ.modelR.s2, cen=0, by=0.2)
typ.pred.rain.s3 <- crosspred(typ.cb.rain.s3, typ.modelR.s3, cen=0, by=0.2)
typ.pred.rain.s4 <- crosspred(typ.cb.rain.s4, typ.modelR.s4, cen=0, by=0.2)
typ.pred.rain.s5 <- crosspred(typ.cb.rain.s5, typ.modelR.s5, cen=0, by=0.2)
typ.pred.rain.s6 <- crosspred(typ.cb.rain.s6, typ.modelR.s6, cen=0, by=0.2)

typ.pred.temp.s1 <- crosspred(typ.cb.temp.s1, typ.modelT.s1, cen=23, by=0.2)
typ.pred.temp.s2 <- crosspred(typ.cb.temp.s2, typ.modelT.s2, cen=23, by=0.2)
typ.pred.temp.s3 <- crosspred(typ.cb.temp.s3, typ.modelT.s3, cen=23, by=0.2)
typ.pred.temp.s4 <- crosspred(typ.cb.temp.s4, typ.modelT.s4, cen=23, by=0.2)
typ.pred.temp.s5 <- crosspred(typ.cb.temp.s5, typ.modelT.s5, cen=23, by=0.2)
typ.pred.temp.s6 <- crosspred(typ.cb.temp.s6, typ.modelT.s6, cen=23, by=0.2)
    
par(mfrow=c(2,6))
plot(typ.pred.rain.s1,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="A")
plot(typ.pred.rain.s2,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="C")
plot(typ.pred.rain.s3,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="E")
plot(typ.pred.rain.s4,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="G")
plot(typ.pred.rain.s5,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="I")
plot(typ.pred.rain.s6,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="K")
    
plot(typ.pred.temp.s1,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR oftyphoid",cex.lab=1.1,cex.axis=1.1,main="B")
plot(typ.pred.temp.s2,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="D")
plot(typ.pred.temp.s3,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="F")
plot(typ.pred.temp.s4,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="H")
plot(typ.pred.temp.s5,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="J")
plot(typ.pred.temp.s6,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="L")

Descriptive stats of the study population
mean(mo.dlnmN$incid_seaX) 
## [1] 6.795455
ci(mo.dlnmN$incid_seaX)
##   Estimate   CI lower   CI upper Std. Error 
##  6.7954545  6.1608513  7.4300578  0.3207921
mean(mo.dlnmT$incid_seaX)
## [1] 3.633333
ci(mo.dlnmT$incid_seaX)
##   Estimate   CI lower   CI upper Std. Error 
##  3.6333333  2.9011782  4.3654884  0.3658954